function result = canonical_correlation_analysis(data,freq,nHarmonics,stimTime,Fs)
%---input---
%data:数据，size1通道，size2数据点
%freq：所有刺激频率，单位Hz
%nHarmonics：谐波次数
%stimTime：窗长，单位s
%Freq：采样率，单位Hz
%---output---
%result：识别结果
data = data';
N = round(stimTime*Fs);
t = (1:N)/Fs;

%% reference signals for CCA
for f = 1:length(freq)
    for nh = 1:nHarmonics
        s1 = sin(2*pi*nh*freq(f)*t);
        s2 = cos(2*pi*nh*freq(f)*t);
        refData(:,2*nh-1:2*nh,f) = cat(2,s1',s2');
    end
end

%% CCA Algorithm
bpData = data(:,:);
cc=[0,0,0,0];
Kcc=[0,0,0,0];
KKcc=[1,1,1,1];
parfor j=1:length(freq)
    [~,~,D] = canoncorr(bpData, refData(:,:,j));
    cc(j)=sum(D)/length(D);
end
cc=(cc-Kcc).*KKcc;
% target identification
result = find(cc==max(cc));
end

